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Abstract 

The  scattering  of  electromagnetic  waves  from  complex  object  coated  with  lossy  materials  with 
negative  and  positive  permitivity  embedded  in  a  layered  medium  is  analysed  via  a  new  formulation, 
and  a  new  set  of  algorithms  to  implement  these  formulations  is  introduced.  This  new  technology  is 
applied  here  to  phase  shifted  photomasks.  The  results  for  various  wavelengths  and  mask 
thicknesses  have  been  considered  carefully  and  are  reported  here. 


1 .  Introduction 

The  problem  of  accurately  computing  the  electromagnetic  scattering  from  lossy,  dispersive  objects 
of  complex  shapes  embedded  in  a  layered  dispersive  medium  has  been  outstanding  for  quite  some 
time,  despite  the  large  range  of  potential  military  applications.  These  applications  include  the 
detection  and  identification  of  targets  located  under  a  foliage  covered  earth  or  airborne  targets  over 
the  horizon.  The  difficulties  associated  with  such  problems  in  part  originate  from  the  need  to  model 
the  shapes  and  inhomogeneities  of  the  objects  accurately  while  simultaneously  keeping 
computational  complexity  to  a  manageable  level. 

The  classic  finite-difference  time-domain  (FDTD)  method  and  finite-element  method  (FEM)  are 
inadequate  for  addressing  the  above  mentioned  difficulties.  This  is  because  FDTD  requires  a 
regular  computational  grid  and  is  therefore  not  suitable  for  modeling  curved  surfaces  accurately.  On 
the  other  hand,  FEM  can  model  curved  surfaces  accurately,  but  requires  the  inversion  of  a  large 
matrix,  which  is  computationally  expensive  for  large  problems. 

For  mask  and  wafer  analysis,  in  part  due  to  the  periodicity  of  some  structures  like  contact  holes  and 
in  part  due  to  features  complexity,  we  have  elected  to  employ  a  new  formulation  of  the  basic 
problem  to  enhance  the  necessary  hybrid  method  for  such  complex  problems.  The  “single  integral 
equation”  method  was  originally  based  on  the  scalar  Green’s  function.  We  have  reformulated  the 
“single  integral  equation”  method  based  on  the  dyadic  Green’s  functions  for  the  layered  wafer 
substrate.  This  approach  allows  each  feature  on  the  wafer  to  be  treated  as  a  localized  object 
embedded  in  a  layered  medium,  resulting  in  a  smaller  total  number  of  unknowns  along  the  profile 
(or  surface)  of  the  embedded  scattering  object  only.  Furthermore,  the  dyadic  Green’s  functions  for  a 


given  layered  structure,  wavelength  and  periodicity  need  to  be  computed  only  once.  The  results  are 
stored  in  interpolation  tables  for  use  on  different  feature  shapes.  This  leads  to  a  particularly  efficient 
methodology  for  investigating  the  effect  of  change  of  the  feature  shape  on  the  diffraction  efficiency 
for  fixed  layered  structure,  wavelength  and  periodicity,  presenting  a  very  powerful  tool  for 
generating  libraries  for  the  desired  analysis. 

The  embedded  object  representing  the  photoresist  feature  under  study  behaves  like  a  cavity 
trapping  an  oscillatory  wave,  enabling  the  generation  of  a  resonance.  This  resonance  expresses 
itself  in  various  ways,  such  as  a  vanishing  eigenvalue  or  a  discontinuous  derivative,  leading  to  a 
potentially  non-unique  solution  for  some  specific  frequencies.  The  “single  integral  equation”  can  be 
formulated  in  terms  of  either  the  electric  field  or  the  magnetic  field  expression,  the  two 
formulations  being  completely  equivalent  in  determining  the  fields.  However,  a  unique  solution  of 
all  frequencies,  including  the  resonant  frequencies  of  the  embedded  object  representing  the  wafer 
feature  under  study,  is  obtained  by  employing  a  linear  combination  of  the  magnetic-field  and 
electric-field  expressions. 

We  have  demonstrated  the  formulation  stability  for  all  frequencies,  as  well  as  the  potential  pitfalls 
associated  with  naive  coding  ignoring  this  most  subtle  point.  We  demonstrate  these  considerations 
on  several  test  pattern  of  photoresist  structure. 

This  report  is  organized  as  follows: 

1 .  In  the  remainder  of  this  Section,  a  program  flowchart  and  vertical  versus  horizontal  critical 
dimension  results  are  presented  (numerical  date,  intensity  plots,  and  aerial  images). 

2.  In  Section  2,  a  basic  formulation  of  the  problem  is  presented. 

3.  In  the  Appendix,  a  preprint  of  the  paper,  “Application  of  the  Hybrid  Finite-Difference 
Time-Domain  Method  to  Modeling  Curved  Surfaces  in  Three-Dimensional  Lithography 
Simulation”  is  presented  for  background  on  the  present  approach. 
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PROGRAM  FLOW  CHART 
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2.  Basic  Formulation 


2. 1 .  Scope  of  Work 

The  work  involves  the  numerical  simulation  of  the  effects  of  polarization  on  aerial-image  quality  in 
157-nm  lithography.  The  finite-difference  time-domain  (FDTD)  algorithm  is  being  used  to  compute 
the  electromagnetic  field  distribution  in  the  near  zone  of  a  chromium  photomask  illuminated  by  a 
plane  wave  with  arb  itrary  polarization  and  angle  of  incidence.  The  corresponding  aerial  image  in 
the  far  zone  is  then  being  computed.  The  effects  of  partial  coherence  is  being  taken  into  account  by 
incoherent  superposition  of  the  aerial  images  corresponding  to  different  angles  of  incidence. 


2.2.  Innovative  Claim 

Three  innovative  features  distinguish  our  approach  from  other  approaches: 

1 .  oBiique  incidence  is  handled  by  the  use  of  a  Huygens  surface  excitation,  in  which  the  exact 
numerical  dispersion  relations  and  numerical  reflection  coefficients  for  the  finite-difference 
mesh  are  employed  to  compute  the  steady-state  field  distribution  within  the  multilayered 
photomask  blank,  for  each  Fourier  component  of  the  incident  pulse. 

2.  A  second-order  accurate  finite-difference  updating  equation  for  the  electric  field  is  used  for 
the  chromium  region  of  the  photomask  structure,  which  is  modeled  by  a  plasma  permittivity 
function. 

3.  A  near-to-far  field  transformation  based  on  multilayered  media  Green's  functions  is  used  to 
compute  the  aerial  image  in  the  far  zone  from  the  field  distribution  in  the  near  zone,  without 
assuming  periodic  boundary  conditions. 

The  mathematical  formulation  of  these  innovative  features  is  discussed  below. 


2.3.  Mathematical  Formulation 

In  this  section,  we  discuss  the  mathematical  formulation  of  the  above  innovative  features  of  the 
software  that  is  currently  being  developed  for  the  simulation  work. 


2.3.1  Huygens  Surface  Excitation 

In  most  existing  FDTD  software  packages  for  lithography  simulation,  the  excitation  is  applied  to 
the  uppermost  surface  of  the  computational  domain  only.  Such  an  approach  is  strictly  applicable 
only  to  normal  incidence,  since,  in  the  case  of  oblique  incidence,  the  excitation  enters  the 
computational  domain  from  the  sides  also.  In  the  general  case,  it  is  necessary  to  apply  the  excitation 
to  all  six  sides  of  the  computational  domain,  or,  equivalently,  to  a  Huygens  surface  completely 
enclosing  the  mask  feature. 
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The  technique  of  Huygens  surface  excitation  was  developed  by  Holland  and  Williams  [1]  for 
objects  embedded  in  empty  space.  For  lithography  simulation,  however,  it  is  necessary  to  take  into 
account  the  presence  of  the  substrate,  which  in  our  case  is  the  multilayered  photomask  blank 
extending  essentially  to  infinity  horizontally.  Here,  it  is  necessary  to  generalize  Holland's  method  to 
take  into  account  the  multiple  reflections  of  the  incident  plane  wave  inside  the  multilayered 
substrate. 

Although  the  propagation  of  a  plane  wave  in  a  multilayered  medium  is  a  well-known  problem  in 
classical  optics,  one  cannot  apply  the  weli-known  results  of  thin-film  optics  directly  to  FDTD 
computation.  The  reason  is  that,  in  the  finite-difference  method,  one  deals  with  a  discrete  lattice 
space,  rather  than  a  continuous  space  on  which  the  thin-film  optics  solution  is  based.  Indeed,  a 
plane  wave  traveling  in  a  lattice  space  obeys  a  different  dispersion  relation,  that  is,  has  a  different 
phase  velocity,  than  one  traveling  in  continuous  space.  Furthermore,  the  Fresnel  coefficients  for  a 
lattice  space  are  different  from  those  for  a  continuous  space.  In  order  to  generalize  Holland's 
method  to  lithography  simulation,  one  must  first  determine  the  correct  dispersion  relations  and 
Fresnel  coefficients  for  a  multilayered  lattice  space.  Furthermore,  one  must  consider  a  plasma 
lattice  space,  since  the  materials  used  in  157-nm  lithography  can  have  negative  real  permittivity, 
that  is,  a  refractive  index  n  =  Ur  +  Jk  for  which  ac  >  Ur,  which  requires  the  use  of  a  plasma 
permittivity  function. 


2.3. 1 . 1  Dispersion  Relation  for  a  Plasma  Lattice  Space 


The  dispersion  relation  for  a  plane  wave  traveling  in  an  ordinary  (non-plasma)  lattice  space  is  well 
known.  It  is  included  here  for  reference: 

sin2(^)  sin2(^^)  ^  sin^ 

(Ax)2  (Ay)2  ^  {Azif  ~  (Af)2  ’  ^  ’ 

where  e\  >  0  is  the  permittivity  of  the  ordinary  dielectric  medium. 


The  dispersion  relation  for  a  plasma  lattice  space,  however,  is  not  so  well  known.  We  have  derived 
the  following  new  result  for  a  plasma  lattice  space  with  plasma  parameters  ujp  and  Oc' 
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The  dispersion  relations  Eqs.  (1)  and  (2)  are  used  to  determine  the  Fresnel  coefficients  for  the 
interface  between  two  dissimilar  lattice  spaces,  as  discussed  in  the  next  section. 


2.3. 1 .2  Fresnel  Coeehcients  for  a  Plasma  Lattice  Space 


Consider  the  interface  between  two  dissimilar  lattice  spaces  as  shown  in  Fig.  1 ,  where  the  region 
2  >  0  is  an  ordinary  lattice  space  with  permittivity  ei  and  the  region  2  <  0  is  a  plasma  lattice 
space  with  parameters  ujp  and  Uc.  Let  a  plane  wave  be  incident  obliquely  from  the  upper  region. 
For  simplicity,  we  consider  the  two-dimensional  case,  where  the  wavevector  of  the  incident  wave  is 
(/Cy,  kz)  with  =  0.  In  order  to  determine  the  Fresnel  coefficients  R  and  T  for  the  interface,  we 
assume  the  following  forms  for  the  electric  and  magnetic  fields  in  the  upper  and  lower  lattice 
spaces,  respectively: 
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where  the  coefficients  Ai  and  Bi  are  determined  by  the  finite-difference  field  equations  for  the 
upper  {i  =  1)  and  lower  {i  =  2)  lattice  spaces,  respectively.  Also,  kzi  and  kz2  are  determined  by 
the  dispersion  relations  Eqs  (1)  and  (2),  respectively.  From  Fig.  1,  it  can  be  seen  that  there  is  an 
electric-field  node  right  at  the  interface  z  =  0.  The  electric  fie'd  ^  at  this  node  is  given  by 

neither  Eq.  (5)  nor  Eq.  (6).  Instead,  we  assume  it  to  have  the  following  form 


pn  _  p  j{ujnAt-\-pkyAy) 

0  - 


(7) 


where  F  is  another  unknown. 

Altogether,  there  are  three  unknowns  in  the  above  problem,  namely,  R,  T  and  F.  To  solve  for 
the.se  unknowns,  one  needs  three  equations.  Two  of  these  equations  are  the  finite-difference 
updating  equations 

for  the  magnetic-field  nodes  marked  Ni  and  N2  in  Fig.  1,  respectively. 
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The  third  equation  is  obtained  by  apply  Ampere's  Law  to  the  circuit  encircling  the  electric-field 
node  at  the  interface  z  =  0  and  bounded  by  the  magnetic-field  nodes  Ni  and  N2'. 
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Eqs.  (8)  t9„(10)  allow  the  Fresnel  coefficients  R  and  T,  as  well  as  the  interfacial  unknown  F,  to  be 
solved.  The  numerical  results  are  illustrated  in  Fig.  2,  which  shows  the  computed  coefficients  for  a 
uniform  lattice  spacing  of  Ay  =  Azi  =  Az2  =  A/15,  which  is  the  typical  mesh  spacing  used  in 
an  FDTD  simulation.  Also  shown  in  Fig.  2  are  the  Fresnel  coefficients  given  by  thin-film  optics  for 
continuous  space.  It  can  be  seen  that,  at  a  mesh  spacing  of  A/15,  there  are  significant  differences 
between  the  coefficients  for  the  lattice  space  and  and  those  for  continuous  space.  To  verify  the 
correctness  of  the  above  formulation,  we  repeated  the  calculation  for  a  uniform  mesh  spacing  of 
A/100.  For  such  a  small  mesh  spacing,  we  expect  the  results  for  the  lattice  space  to  converge  to 
those  for  continuous  space.  This  is  indeed  the  case,  as  shown  in  Fig.  3. 


The  Fresnel  coefficients  for  the  lattice  space  must  be  separately  computed  by  the  above  procedure 
for  each  of  the  several  dozen  Fourier  components  needed  to  characterize  a  finite-duration,  pulse 
excitation.  The  corresponding  multiply  reflected  waves  of  different  frequencies  inside  the  lattice 
substrate  must  afterwards  be  re-assembled  to  produce  the  correct  time-domain  excitation  waveform 
at  each  point  on  the  Huygens  surface. 


2.3.2  Second-Order  Accurate  Discretization  for  the  Plasma  Region 


The  plasma  model  of  a  dispersive  medium  was  introduced  by  Luebbers,  et  al.  [2].  However,  the 
updating  equation  for  the  electric  field  given  by  these  authors  is  only  first-order  accurate  in  time, 
due  in  part  to  the  use  of  the  rectangular  rule  for  the  convolution  integral.  In  our  implementation, 
however,  the  electric  field  is  assumed  to  vary  linearly  with  time  between  successive  time  steps  and 
the  resulting  convolution  integral  is  integrated  exactly.  This  way,  the  following  second-order 
accurate  updating  equation  for  the  electric  field  has  been  derived  [3]: 
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where  Gq  and  Cti  are  given  by  Eqs.  (3)  and  (4),  respectively,  and  the  vector  is  defined  as 
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which  can  be  computed  recursively  by 

[ai  (l  +  E"-^  +  .  (13) 

Numerical  experiments  have  indicated  that  Eq.  (11)  is  more  accurate  than  the  first-order  accurate 
updating  equation  given  by  Luebbers,  et  al.  [2]  which  is  used  in  most  existing  FDTD  software 
packages. 

2.3.3  Near-to-Far  Field  Transformation 

An  FDTI^..pomputation  only  gives  the  fields  in  the  near  zone  of  the  mask  feature.  To  obtain  the 
aerial  image  in  the  far  zone,  one  must  extrapolate  the  near-zone  fields  to  the  far  zone  using  the 
Kirchhoff  integral  representation.  In  most  existing  FDTD  software  packages,  this  is  accomplished 
by  integrating  the  fields  over  a  horizontal  surface  lying  just  outside  the  photomask  blank  on  the 
projection-lens  side.  However,  since  the  photomask  blank  extends  essentially  to  infinity 
horizontally,  it  would  in  principle  be  necessary  to  integrate  over  an  infinite  horizontal  surface  to 
compute  the  far-zone  fields,  which  is  clearly  impractical  unless  the  geometry  is  artificially 
truncated  by  assuming  periodic  boundary  conditions,  as  is  usually  done. 

We  have  recently  developed  a  near-to-far  field  transformation  which  does  not  assume  periodic 
boundary  conditions  but  yet  requires  only  integration  over  a  finite  surface  S,  which  we  shall  call 
the  Kirchhoff  surface,  enclosing  the  mask  feature  of  interest  [4].  In  this  method,  the 
electromagnetic  vector  potentials  are  expressed  in  terms  of  multilayered  media  Green's  functions 
for  the  photomask  substrate, 

A^^^(r)  =  //o  /  £  (r|r')  ■  Je(r ')  ^5"  ,  (14) 

=  eo^[di(r|r')-Jm(r')  +  VCM(r|r')Jn,..(r')]  ,  (15) 

where  Qe  and  Gm  are  multilayered  media  Green's  functions, 

&(r|r')  =  (xx  +  yy)Gf,(r|r')  +  zxGf,(r|r')  +  zyGf,(r|r')  +  zzG£(r|r')  ,  (16) 
el,  (r|r')  =  (xx  +  yy)G"(r|r')  +  zxG"(r|r')  +  zyG"(r|r')  +  zzG"(r|r')  .  (17) 

Using  the  technique  of  asymptotic  evaluation  of  integrals,  we  have  derived  closed-form  expression 
for  the  multilayered  media  Green's  functions  for  the  far-zone  fields: 
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and  di  is  the  thickness  of  the  chromium  layer. 
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2.4.  Algorithmic  Implementation 

In  this  section,  we  discuss  how  the  mathematical  formulation  discussed  in  Section  3  has  been 
implemented  in  computer  code. 

For  a  given  angle  of  incidence  and  polarization,  the  program  reads  in  the  geometry  description 
contained  in  an  input  file  and  then  proceeds  to  set  up  a  finite-difference  computational  mesh  and 
calculate  the  various  parameters  used  in  the  updating  equations  for  the  electric  and  magnetic  fields 
at  each  node.  Next,  the  incident  pulse  waveform  is  decomposed  into  Fourier  components.  For  each 
Fourier  component,  the  dispersion  relations  in  the  various  material  regions  are  evaluated  using  Eqs. 
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(1)  and  (2).  The  results  are  substituted  into  the  system  of  equations  Eqs.  (8)  to  (10)  for  each 
material  interface  in  the  multilayered  photomask  substrate.  These  equations  are  solved  for  the 
Fresnel  coefficients  R  and  T  associated  with  that  material  interface.  Using  these  coefficients,  the 
steady-state  field  distribution  in  the  photomask  substrate  for  the  particular  Fourier  component  is 
computed  and  stored  at  all  the  points  on  the  Huygens  surface.  This  process  is  repeated  for  all  the 
Fourier  components  of  the  incident  waveform.  Afterwards,  the  stored  frequency  reponse  at  each 
point  on  the  Huygens  surface  is  converted  into  a  time-domain  waveform  by  inverse  Fourier 
transformation.  This  completes  the  initialization  of  the  incident  excitation  on  the  Huygens  surface. 

The  program  then  enters  the  main  FDTD  time-marching  loop.  The  fields  in  the  entire 
computational  domain  are  initialized  to  zero.  At  the  start  of  the  n“'  time  step,  the  magnetic  field 
components  at  all  the  magnetic-field  nodes  in  the  domain  are  updated  to  their  values  at  time  step 
n  -f  The  incident  excitation  for  time  step  rt  -h  |  is  then  applied  to  the  magnetic-field  nodes  on 
the  Huygens  surface.  Next,  the  electric  field  components  at  all  the  electric-field  nodes  in  the 
domain  are  updated  to  their  values  at  time  step  ti  +  1.  For  electric-field  nodes  lying  in  the 
chromium -region  of  the  photomask  structure,  the  second-order  accurate  equation  Eq.  (11)  is  used, 
while  for  the  remaining  electric-field  nodes,  the  standard  FDTD  updating  equation  is  used.  The 
incident  excitation  for  time  step  +  1  is  then  applied  to  the  electric-field  nodes  on  the  Huygens 
surface.  The  nodes  on  the  outermost  boundaries  of  the  computational  domain  have  yet  to  be 
updated.  An  absorbing  boundary  condition  (ABC)  is  needed  for  this  purpose.  Various  ABCs  are 
available  in  the  literature,  from  simple  but  less  accurate  ones,  such  as  the  Mur  ABC,  to 
sophisticated  and  highly  accurate  ones,  such  as  Berenger's  Perfectly  Matched  Layer.  In  our  code, 
we  have  used  the  first-order  Higdon  ABC,  which  has  been  shown  to  give  accurate  results  provided 
that  the  mask  feature  under  investigation  is  separated  from  the  outermost  boundaries  by  a  distance 
of  at  least  one  wavelength.  After  updating  the  fields  on  the  outermost  boundaries  using  the  Higdon 
ABC,  the  FDTD  loop  for  the  time  step  is  completed.  Before  proceeding  with  the  next  time  step, 
however,  the  discrete  Fourier  transform  of  the  field  at  each  node  on  the  Kirchhoff  surface  S  is 
computed  recursively  and  saved.  This  way,  upon  completion  of  the  FDTD  loop  for  all  the  time 
steps,  we  shall  immediately  have  available  the  steady-state  response  of  the  system  at  each  point  on 
the  Kirchhoff  surface,  which  has  been  used  in  the  subsequent  near-to-far  field  transformation. 

Since  a  pulse  of  finite-duration  is  used  for  excitation,  the  fields  in  the  computational  domain 
eventually  decays  with  time.  When  the  fields  have  decayed  to  a  sufficiently  low  level,  the  FDTD 
time-marching  loop  is  exited.  The  computed  steady-state  fields  on  the  Kirchhoff  surface  are  then 
written  to  an  output  file  and  the  FDTD  program  terminates.  Another  program  is  executed  to 
perform  the  near-to-far  field  transformation  on  the  output  of  the  FDTD  program.  The  program 
computes  the  far-zone  fields  at  selected  points  on  the  entrance  pupil  of  the  projection  lens  using 
Eqs.  (14)  to  (29).  Then,  the  fields  on  the  entrance  pupil  are  propagated  to  the  exit  pupil,  and  from 
there  to  the  wafer  plane,  by  using  the  vector  aerial  image  model  [5]. 

The  above  calculations  is  being  repeated  for  different  angles  of  incidence  in  order  to  simulate  the 
effects  of  partial  coherence.  This  process  has  been  automated  by  the  use  of  shell  scripts  to  control 
program  execution.  It  should  also  be  pointed  out  that  our  data  structure  has  been  designed  in  such  a 
way  as  to  facilitate  future  upgrading  of  the  FDTD  algorithm  used  for  this  work  to  the  more 
powerful  hybrid  finite-element  and  finite-difference  time-domain  algorithm  [6]. 
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Fig.  I .  Interface  between  an  ordinary  dielectric  lattice  space  and  a  plasma  lattice 
space. 
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Reflection  coefficient  (R)  Reflection  coefficient  (R) 


I'ig.  2.  Reflection  and  transmission  coefficients  for  the  interface  between  air 
and  a  plasma  material  with  n  =  0.85  +  2.01/.  Dashed  lines  arc  for 
discrete  lattice  spaces  with  a  unifoiTn  lattice  spacing  of  A  “  ^15.  Solid 
lines  arc  for  continuous  space 


Mg,  3.  Same  as  Fig.  2,  but  for  a  uniform  lattice  spacing  of  A  ”  ^1 00. 
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Appendix 


Application  of  the  Hybrid  Finite-Difference 
Time-Domain  Method  to 
Modeling  Curved  Surfaces  in 
Three-Dimensional  Lithography  Simulation 


1.  Introduction 

Electromagnetic  scattering  from  nonplanar  topography  on  the  mask  or  wafer  can  have  a  significant 
impact  on  photomask  performance,  linewidth  control  or  alignment  precision.  Computer  simulation 
is  a  cost-effective  way  to  assess  the  importance  of  such  electromagnetic  scattering  effects  in 
photolithography.  With  the  scaling  of  devices  to  smaller  dimensions,  greater  demands  are  placed 
on  the  accuracy  of  the  mathematical  models  used  in  the  electromagnetic  simulators  and  the 
efficiency  of  their  numerical  implementation. 

Over  the  past  many  years,  the  finite-difference  time-domain  method  (FDTD)  has  become  the 
prevalent  method  for  solving  electromagnetic  scattering  problems.  Not  only  is  it  generally 
applicable  to  arbitrary  geometries,  but  also  it  is  the  most  efficient  algorithm  available,  by  yielding 
useful  field  information  at  Aspace  points  and  m  time  points  in  a  total  of  only  0{m.N)  operations. 
Thus,  for  example,  a  single  time-domain  simulation  of  the  response  of  an  electromagnetic  system 
to  a  finite-duration  incident  pulse  can  yield  the  steady-state  results  for  a  large  number  of  different 
frequencies  by  Fourier  transformation. 

The  main  disadvantage  of  FDTD  is  its  inefficiency  in  modeling  curved  surfaces  accurately,  since 
the  regular  finite-difference  mesh  used  in  FDTD  requires  that  curved  surfaces  be  approximated  by 
staircase  models.' '^o  achieve  accurate  results  using  the  staircase  model,  one  usually  has  to  use  a 
very  fine  mesh  and,  therefore,  also  a  very  small  time  step  due  to  the  stability  criterion. 

Recently,  a  hybrid-FDTD  method  appeared  in  the  literature  [1]  which  combines  the  flexibility  of 
the  finite-element  method  (FEM)  in  modeling  curved  surfaces  accurately  with  the  computational 
efficiency  of  FDTD.  It  appears  that  this  hybrid-FDTD  method,  when  used  in  conjunction  with 
high-performance  absorbing  boundary  conditions,  is  the  ideal  tool  for  solving  three-dimensional 
electromagnetic  scattering  problems  arising  in  lithography  simulation  accurately  and  efficiently. 

However,  the  hybrid-FDTD  method  in  its  original  form  is  not  suitable  for  DUV  lithography 
simulation.  This  is  because  the  original  formulation  cannot  handle  lossy  materials,  especially  those 
with  negative  dielectric  constants,  such  as  chromium  and  silicon,  which  are  commonplace  in  DUV 
lithography.  The  goal  of  this  paper  is  to  extend  the  orignal  hybrid-FDTD  formulation  to  handle 
lossy  materials. 

After  reviewing  the  original  hybrid-FDTD  formulation  in  Section  2.1,  we  discuss  its  extensions  to 
lossy  materials  with  positive  and  negative  dielectric  constants  separately  in  Sections  2.2  and  2.3. 
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The  issues  of  mesh  generation  are  discussed  briefly  in  Section  3.1.  This  is  followed  in  Section  3.2 
by  a  presentation  of  numerical  results  comparing  the  hybrid-FDTD  method  and  standard  FDTD  for 
various  dielectric  materials. 


2.  The  Hybrid-FDTD  Method 
2. 1  The  Original  Method 

The  hybrid-FDTD  method  was  originally  developed  by  Wu  and  Itoh  [1]  for  lossless  dielectric 
objects.  Consider  a  dielectric  object  bounded  by  a  curved  surface  S  as  shown  in  Fig.  1  The 
computational  domain  is  divided  into  two  overlapping  regions:  (i)  A  regular  finite-difference 
region  Hi  spanning  the  interior  and  exterior  of  the  object  at  some  distance  from  the  surface  S,  and 
(ii)  an  irregular  finite-element  region  H2  spanning  the  immediate  vicinity  of  S  on  both  of  its  sides. 
The  two  regions  overlap  in  a  single  layer  of  finite-difference  cells  bounded  by  staircase  surfaces  F 1 
and  r2on'each  side  of  S. 

Suppose  the  electric  field  E”is  known  everywhere  at  time  step  n.  Using  the  standard  Yee 
algorithm  [2],  the  magnetic  field  at  time  step  n  -h  ^  in  the  regular  region  Hi,  including  the 

overlap  region,  can  be  updated.  This  in  turn  allows  the  electric  field  in  Hi,  including  the 

boundary  Fi,  to  be  updated.  The  electric  field  E""^^  in  the  irregular  region  H2  is  then  updated  by 
solving  the  weighted-residual  problem 

— /  eE“-EdH  =  -  /  -VxE“-VxEdH  ,  (1) 

of/  .la  .In  // 

using  the  electric  field  E"'^^on  Fi  as  boundary  condition  and  the  previous  electric  fields  E"and 
E"”^  in  H2  as  initial  condition.  Once  the  electric  field  at  time  step  n -|- 1  becomes  available 
everywhere,  the  time-marching  can  be  continued  for  the  next  time  step. 

To  solve  the  weighted-residual  problem  Eq.  (1),  the  irregular  region  H2  is  subdivided  into  many 
small  tetrahedral  elements  and  the  electric  field  in  the  elements  is  expanded  in  Whitney  vector  basis 
functions  W,  [3],  with  the  electric  field  components  Ej  along  the  element  edges  jas  the  expansion 
coefficients.  Next,  the  time  derivative  in  Eq.  (1)  is  approximated  by  the  central  difference  operator 
and  the  Newmark-Beta  method  [1]  is  applied  to  obtain  an  unconditionally  stable,  second-order 
accurate,  implicit  time-marching  scheme, 

Acl  +  £"«  =  2  Ac|  -  E”  -  Ac]  +  E”-'  .  (2) 

Here  the  matrices  [C]  and  [D]  are  given  by 
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ich  = 

/  eoCr-W,  •  W,  dQ  , 

./n 

(3) 

= 

[  -V  xWi-V  xW.dQ  , 

Jq 

(4) 

where  Cr  is  the  dielectric  constant  and  is  the  vector  of  expansion  coefficients  at  time  step  n. 

The  above  formulation  of  Wu  and  Itoh  is  applicable  only  to  lossless  dielectric  objects  with  positive 
dielectric  constants  tr.  The  objects  encountered  in  lithography  simulation,  however,  are  often  lossy 
and,  furthermore,  may  have  negative  dielectric  constants,  especially  at  DUV  wavelengths.  We  next 
discuss  the  extensions  of  the  formulation  of  Wu  and  Itoh  to  lossy  media  with  positive  and  negative 
dielectric  constants  separately. 

2.2  Extension  to  Lossy  Medium  with  Positive  Dielectric  Constant 

A  material-,- with  positive  dielectric  constant  is  one  for  which  the  real  part  n  of  its  complex 
refractive  index  is  greater  than  the  imaginary  part  k.  Such  is  the  case  for  weakly  absorbing 
materials  such  as  photoresist  and  silicon  nitride.  This  kind  of  material  can  be  modeled  in  the  time 
domain  by  adding  a  conductivity  term  to  the  electric-field  updating  equation.  The  dielectric 
constant  Cr  and  conductivity  cr  are  related  to  the  complex  refractive  index  n  -f  Jk  by 

Cr  =  ,  (5) 

a  -  2nK.iJoeQer  ,  (6) 


where  loq  is  the  frequency  of  interest. 

By  using  exponential  time  stepping  [4],  the  updating  equation  for  the  electric  field  in  the  irregular 
region  becomes 

(|C+|  +  =  2  (iCol  -  B’"'  ,  (7) 

where  the  matrices  [C±]  and  [Co]  are  given  by 

[C±]rj  =  ^  eoe.  exp  (±^)  W,;  •  W,-  dQ  ,  (8) 

[Co]i,  =  ^  eoer  cosh  W,;  •  W,-  dQ  .  (9) 

The  time-marching  scheme  Eq.  (7)  is  unconditionally  stable  as  long  as  Cr  >  0. 

2.3  Extension  to  Lossy  Medium  with  Negative  Dielectric  Constant 

A  material  for  which  n  <  k  has  a  negative  dielectric  constant  according  to  Eq.  (5).  Silicon, 
chromium  and  tungsten  are  examples  of  such  materials  at  DUV  wavelengths.  This  type  of  material 


65 


can  be  modeled  in  the  time  domain  by  an  unmagnetized  plasma  [5],  with  a  complex  dielectric 
function 


e^(w) 


1  - 


uj{uj  +  jiyc) 


(10) 


The  plasma  frequency  LOp  and  collision  frequency  I'c  appearing  in  Eq.  (10)  are  related  to  the 
complex  refractive  index  of  the  material  at  the  frequency  loq  of  interest  by 


U!, 


P  _ 


Wo 

UlQ 


1  +  k"  - 
2nK 


n^  + 


An  K 


2  ,.2 


1  + 


1  +  —  n? 


(11) 

(12) 


The  weighted-residual  problem  to  be  solved  in  the  irregular  plasma  region  is 


/  eoE“  •  E  dfi  =  -  /  -  V  X  E“  •  V  X  E  dfi  -  /  ■  E  dQ 

ot^  -In  .In  //  Jn  ‘ 

+  ['  dt!  f  •  E{t  -  t!)  dQ  ,  (13) 

.1—00  Jo. 

in  which  a  convolution  term  appears  due  to  the  frequency  dependence  of  the  plasma  dielectric 
function  Eq.  (lO)Expanding  the  electric  field  E  in  Whitney  basis  functions  and  applying  the 
Newmark-Beta  method,  we  obtain  the  following  second-order  accurate,  implicit  time-marching 
scheme: 


-  (f^i + ^  (I®) + h"’]})  *  ‘  •  (1-*) 


where  [C]  is  given  by  Eq.  (3)  with  =  1 


and  the  matrix 


is  given  by 


) 


Wi  ■  Wj  dD.  , 


The  vector  appearing  in  Eq.  (14)  is  defined  as 


7n=l 


(15) 


(16) 


where  the  matrix 


is  given  by 
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A  777 


1  =  I- 

./n 


cosh(i/cAi)  —  1 


UcAt 


Wi  •  W,'  dO. 

V 


(17) 


Using  the  technique  of  Luebbers  et  al.  [5],  the  vector  can  be  computed  recursively  by 


(2)1  ^”-1  _j_  g-t'cAf^n-1 


Xl 


(18) 


The  same  plasma  model  Eq.  (10)  can  be  applied  to  the  regular  FDTD  region  in  the  interior  of  the 
object.  However,  the  traditional  FDTD  implementation  of  this  model  [5]  is  only  first-order  accurate 
in  time,  due  in  part  to  the  use  of  the  rectangular  rule  for  the  convolution  integral.  Instead,  by 
assuming  the  electric  field  to  vary  linearly  with  time  between  successive  time  steps  and  performing 
the  resulting  convolution  integral  exactly,  the  following  second-order  accurate  updating  equation 
for  the  electric  field  in  the  interior  FDTD  region  is  obtained; 
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1  -I-  ^u^At^ao 


I  [l  -  ^ulAt^  (ao  +  E"  -  ^  curl  H”+?](19) 


where 


ao 


Ol 


V  I'cAt  ) 


2 


cosh(t'cAf)  —  1 

KATp 


(20) 

(21) 


The  vector  appearing  in  Eq.  (19)  is  defined  as 


ai  (l  +  E^ 


(22) 


and  can  be  computed  recursively  by 

[ai  (l  +  E""^  +  .  (23) 

To  our  knowledge,  the  second-order  accurate  updating  equation  Eq.  (19)  has  not  appeared  before  in 
the  literature. 


3.  Results 

We  tested  the  formulations  discussed  in  Sections  2.1  to  2.3  with  the  problem  of  electromagnetic 
scattering  by  a  dielectric  sphere,  for  which  the  exact  solution  is  known.  The  diameter  of  the  sphere 
was  0.06  jj,m  and  the  wavelength  was  0.248  fim. 

3.1  Mesh  Generation 
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The  first  step  in  the  computation  was  to  generate  a  high-quality  mesh  for  the  three-dimensional 
irregular  region  Q2  shown  in  Fig.  1 .  After  subdividing  the  surface  of  the  sphere  into  a  large  number 
of  small  triangles,  an  unstructured  tetrahedral  mesh  was  generated  in  Q2  conforming  to  the  surface 
triangulation  of  the  sphere  and  the  staircase  boundaries  Fi  and  r2,  using  our  automatic  mesh- 
generation  software. 

The  quality  Q  of  our  mesh  was  measured  by  the  minimum  of  the  sines  of  all  the  dihedral  angles  in 
the  mesh,  where  0  <  Q  <  1-0.  This  quality  measure  was  chosen  to  bias  against  elements  with  too 
large  ^  18G')  or  too  small  (0min  ~  0°)  dihedral  angles,  which  would  lead  to  poor  accuracy  of 
the  finite-element  interpolation  or  poor  conditioning  of  the  finite-element  matrix,  respectively  [6]. 
It  was  found  that  our  as-generated  mesh  had  a  quality  Q  of  only  0.073,  indicating  the  presence  of 
poorly  shaped  elements  in  the  mesh. 

To  remove  the  poorly  shaped  elements,  we  performed  mesh  improvement  on  the  as-generated  mesh 
in  two  steps.  In  the  first  step,  the  sub-mesh  belonging  to  each  edge  in  the  mesh,  consisting  of  all  the 
tetrahedra  adjacent  to  that  edge,  was  examined.  The  edge  was  deleted  and  replaced  by  one  or  more 
new  edges  if,  by  doing  so,  the  quality  of  the  sub-mesh  was  improved.  In  the  second  step,  the  cluster 
belonging  to  each  node,  consisting  of  all  the  tetrahedra  adjacent  to  that  node  and  its  nearest 
neighbors,  was  examined.  If  the  quality  of  the  cluster  was  below  a  certain  threshold,  the  nodes  in 
the  cluster  were  moved  to  new  positions  which  maximized  the  quality  of  the  cluster.  It  was  found 
that  after  a  single  pass  through  our  mesh-improvement  routine,  the  quality  of  the  mesh  increased  to 
Q  =  0.335,  which  was  deemed  satisfactory  for  our  computation.  The  final  mesh  is  shown  in  Fig.  2 
and  consists  of  1195  nodes,  7158  edges  and  5485  tetrahedra. 

The  above  finite-element  mesh  was  embedded  in  a  regular  finite-difference  mesh  with  20  x  20  X 
20  cells.  For  simplicity,  the  first-order  Higdon  absorbing  boundary  condition  [8]  was  used  on  all 
six  sides  of  the  computational  domain,  which  measured  0.2  /xm  x  0.2  /xm  X  0.2  fim.  A  Huygens 
surface  [9]  located  two  cells  interior  to  the  outermost  boundaries  was  used  to  excite  the  domain 
with  various  Gaus-ian  pulses. 

3.2  Numerical  Results 
3.2.1  Lossless  Dielectric 

We  first  compared  our  results  with  those  of  Wu  and  Itoh,  who  used  a  lossless  dielectric  sphere  of 
refractive  index  U]  =  3.0  and  a  Gaussian  pulse  with  finite  d.c.  content, 

E",  =  xe"M^"01  ,  (24) 

where  Uq  =  33.  The  computed  results  are  shown  in  Fig.  3.  Fig.  3a  shows  the  total  time-domain 
waveform  at  the  center  of  the  sphere,  while  Fig.  3b  shows  the  scattered  waveform  at  a  point  0.09 
jjm  in  front  of  the  sphere.  Also  shown  in  the  figures  are  the  exact  Mie  solution  [7]  and  the  results 
obtained  with  standard  FDTD  using  the  same  mesh  spacing  as  hybrid-FDTD  in  the  regular  region 
Vti,  namely,  20  x  20  x  20  cells,  or  roughly  1/8  of  the  wavelength  rxi  inside  the  dielectric.  The 
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results  of  Figs.  3a  and  3b  are  in  good  agreement  with  those  of  Wu  and  Itoh  [1].  This  verifies  the 
correctness  of  our  computer  program. 

Next,  the  frequency-domain  scattering  cross  section  was  obtained  by  Fourier  transformation  of  the 
corresponding  time-domain  result.  The  hybrid-FDTD  result  is  shown  in  Fig.  3c,  together  with  the 
exact  Mie  solution  and  the  results  obtained  with  standard  FDTD  using  20  x  20  X  20  and  40  X  40  x 
40  cells,  respectively.  It  can  be  seen  that,  compared  with  the  exact  result,  the  hybrid-FDTD  method 
with  a  coarse  mesh  spacing  of  Ai/8  gave  much  better  accuracy  than  FDTD  with  the  same  mesh 
spacing,  and  roughly  the  same  accuracy  as  FDTD  with  the  finer  mesh  spacing  of  Ai/16.  The  small 
discrepancy  between  the  hybrid-FDTD  result  and  the  exact  result  is  due  in  part  to  the  approximate 
absorbing  boundary  conditions  used. 

3.2.2  Lossy  Dielectric  with  Positive  Dielectric  Constant 

We  tested  the  formulation  of  Section  2.2  by  using  a  lossy  dielectric  sphere  of  refractive  index  ni  = 
2.0  +  j0.5-.vThe  time  derivative  of  a  Gaussian  pulse  was  used  for  excitation  to  avoid  introducing  a 
d.c.  offset  into  the  solution. 


E"  = 

me 


where  =  33.  The  computed  time-domain  waveforms  are  shown  in  Figs.  4a  and  4b,  together  with 
the  exact  results  and  the  results  obtained  with  standard  FDTD  using  the  same  mesh  spacing  as 
hybrid-FDTD,  or  roughly  1/12  of  the  wavelength  inside  the  dielectric.  It  can  be  seen  that  the 
hybrid-FDTD  and  FDTD  results  are  both  in  good  agreement  with  the  exact  results,  although  the 
FDTD  results  have  slightly  more  overshooting  at  the  valleys  of  the  waveforms. 


The  result  for  the  Fourier  transformed  scattering  cross  section  is  shown  in  Fig.  4c,  together  with  the 
FDTD  results  obtained  with  a  coarse  and  a  fine  mesh  spacing.  It  can  be  seen  that,  as  in  the  lossless 
dielectric  case,  the  hybrid-FDTD  method  with  a  coarse  mesh  spacing  of  Ai/12  gave  much  better 
accuracy  than  FDTD  with  the  same  mesh  spacing,  and  roughly  the  same  accuracy  as  FDTD  with 
the  finer  mesh  spacing  of  Ai/24. 


3.2.3  Lossy  Dielectric  with  Negative  Dielectric  Constant 


We  tested  the  formulation  of  Section  2.3  by  using  a  lossy  dielectric  sphere  of  refractive  index  n\  = 
0.85  +  j2.01,  which  is  the  refactive  index  of  chromium  at  0.248  /Ltm  [10].  The  time-derivative 
Gaussian  pulse  Eq.  (25)  with  zero  d.c.  content  was  used  for  excitation  to  avoid  the  singularity  of 
the  plasma  dielectric  function  Eq.  (10)  at  w  =  0. 

The  computed  time-domain  waveforms  are  shown  in  Figs.  5a  and  5b,  together  with  the  exact 
results  and  the  results  obtained  with  standard  FDTD  using  the  same  mesh  spacing  as  hybrid-FDTD, 
or  roughly  1/29  of  the  wavelength  inside  the  dielectric  sphere.  It  can  be  seen  that,  whereas  the 
hybrid-FDTD  results  are  in  good  agreement  with  the  exact  results  for  all  times,  the  FDTD  results 
show  marked  departures  from  the  exact  results  at  late  times. 
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The  result  for  the  Fourier  transformed  scattering  cross  section  is  shown  in  Fig.  6,  together  with  the 
FDTD  results  obtained  with  a  coarse,  a  fine  and  a  very  fine  mesh  spacing.  It  can  be  seen  that, 
whereas  the  hybrid-FDTD  result  is  in  good  agreement  with  the  exact  result,  the  FDTD  results  for 
all  three  mesh  spacings,  namely,  Ai/29,  Ai/58  and  Ai/116,  show  large  departures  from  the  exact 
result.  Since  the  second-order  accurate  updating  equation  Eq.  (19)  was  used  for  the  FDTD 
computations,  these  departures  cannot  be  due  to  inaccurate  implementation  of  the  plasma 
dispersion  model  of  Section  2.3,  but,  rather,  must  be  due  to  inaccuracy  of  the  staircase  model  of  the 
spherical  surface  used  in  standard  FDTD.  These  results  highlight  the  need  to  use  the  hybrid-FDTD 
method  to  model  curved  surfaces  accurately  in  the  case  of  lossy  dielectric  materials  with  negative 
dielectric  constants. 


4.  Conclusions 

Extension€--of  the  original  hybrid-FDTD  method  to  handle  lossy  materials  with  positive  and 
negative  dielectric  constants  have  been  discussed  separately.  The  correctness  of  our  computer 
program  has  been  verified  by  comparing  our  computed  results  with  those  in  the  literature  and  with 
the  exact  results.  Our  results  have  shown  that,  for  lossless  dielectric  and  lossy  material  with 
positive  dielectric  constant,  the  hybrid-FDTD  method  is  much  more  accurate  than  standard  FDTD 
when  the  same  mesh  spacing  is  used  in  both  methods,  while  the  two  methods  have  roughly  the 
same  accuracy  when  the  mesh  spacing  used  in  standard  FDTD  is  half  that  used  in  hybrid-FDTD. 
For  lossy  material  with  negative  dielectric  constant,  the  difference  between  the  two  methods  is 
much  more  pronounced.  In  this  case,  the  hybrid-FDTD  method  with  a  mesh  spacing  of  A  is  much 
more  accurate  than  standard  FDTD  even  when  a  mesh  spacing  of  jA  is  used  in  the  latter  method. 
These  results  indicate  that  the  hybrid-FDTD  method  is  far  superior  to  standard  FDTD  for 
lithography  simulation  at  DUV  wavelengths,  where  lossy  materials  with  negative  dielectric 
constants  are  commonplace. 
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Overlap  between  FEM 
and  FDTD  regions 

Huygens  surface 


Fig.  I .  The  hybrid-FDTD  computational  domain  consisting  of  overlapping  FDTD  (Q|)  and 
FEM  (Q2)  regions.  F |  and  r2  are  the  exterior  and  secondary  boundaries,  respectively, 
of  the  overlap  region.  An  incident  wave  is  applied  to  the  Huygens  surface. 
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Fig.  3.  Results  for  a  lossless  dielectric  sf^iere  of  refractive  index  Uj  =  3.0.  (a)  and  (b):  Time 
domain  waveforms  at  center  of  sphere  and  at  a  point  0.09  pm  in  front  of  the  sphere, 
(c):  Radar  cross  section  obtained  by  Fourier  transformation  of  the  time-domain  resul 
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Fig.  4.  Radar  cross  section  of  a  lossy  sphere  of  refractive  index  2.0  -  0.5j,  obtained  by  Fourier 
transformation  of  the  time-domain  results. 
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Fig.  5.  Time-domain  waveforms  for  a  lossy  sphere  of  refractive  index  0.85  -  2.0 ly,  modeled  by 
a  plasma  model,  (a)  Total  wave  at  the  center  of  the  sphere,  (b)  Scattered  wave  at  a  point  0.09  p.m 
in  front  of  the  sphere. 


74 


Fig.  6.  Radar  cross  section  of  a  lossy  sphere  of  refractive  index  0.85  -  2.0  ly,  obtained  by 
Fourier  transformation  of  the  time-domain  results. 
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